Filamentary Diffusion of Cosmic Rays on Small Scales 
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We investigate the diffusion of cosmic rays (CR) close to their sources. Propagating individual 
CRs in purely isotropic turbulent magnetic fields with maximal scale of spatial variations i max , we 
find that CRs diffuse anisotropically at distances r < Z max from their sources. As a result, the CR 
densities around the sources are strongly irregular and show filamentary structures. We determine 
the transition time t» to standard diffusion as t» ~ 10 4 yr (/ max /150pc) /3 (iJ/PeV) '(-Brms/4 uG) 1 , 
with /3 ~ 2 and 7 = 0.25-0.5 for a turbulent field with Kolmogorov power spectrum. We calculate 
the photon emission due to CR interactions with gas and the resulting irregular source images. 



Introduction. — The suggestion that Galactic cosmic 
rays (CR) are accelerated using the energy released in 
supernova (SN) explosions dates back to the 1930's 
This idea was supported initially mainly by the argu- 
ment that SNe inject sufficient energy into the Galaxy to 
maintain the observed CR energy density, while later the 
radio emission observed from SN remnants (SNR) was in- 
terpreted as indication for the acceleration of high-energy 
electrons. Until present, a clear proof for both the accel- 
eration of hadrons and the identity of their sources is still 
missing Alternatively, the required power could be 
provided by acceleration processes which operate at dis- 
tances and time scales larger than for individual SNRs, 
e.g. in superbubbles Q. 

Main obstacle for the identification of CR sources 
is the diffusion of CRs in the Galactic magnetic field 
(GMF), erasing directional information on the position 
of their sources. The GMF has a turbulent compo- 
nent which varies on scales between Z m ; n < 1 AU and 
^max ~ few to 200 pc. Since CRs scatter on inhomo- 
geneities with variation scales comparable to their Lar- 
mor radius, the propagation of Galactic CRs in the GMF 
resembles a random walk and is well described by the dif- 
fusion approximation 0, [B[ . 

However, CRs around young sources do not have time 
to diffuse far away. They should produce an extended 
gamma-ray halo, which can be detected using Cherenkov 
telescopes and gamma-ray satellites. If SNe power the 
Galactic CR population, about ten sources with degree 
extension should be detected at energies E 1 > 100 GeV 
in the Galactic plane assuming the sensitivity of Fermi- 
LAT, while Ref. Q found 18. Most of these sources 
were observed also as extended sources up to energies 



E 1 > lOTcV by the HESS experiment [7fl and have a 
non-spherical shape. Similarly, the Veritas observations 
of Tycho show a clear asymmetric extension of TeV pho- 
tons towards the north of the SNR || . 

The diffusion approximation cannot predict local phe- 
nomena which arise below the CR mean free path A <C 
l m „, when the local configuration of the turbulent field 
must lead to observable imprints Q. Before the present 
study it has remained unclear to which extent the dif- 
fusion approximation is satisfied on intermediate scales 



A <C I < imax, where large-scale fluctuations of the field 
lead to local anisotropics, which in turn can explain the 
irregular images of extended sources found in Refs. 

We study therefore in this letter the diffusion of CRs 
on scales comparable to the coherence length of the tur- 
bulent GMF, O(100 pc). In contrast to earlier studies, we 
calculate the diffusion tensor propagating individual CRs 
in specific realizations of the turbulent magnetic field. 
We find that diffusion is anisotropic even for an isotropic 
random field and can lead to a filamentary structure of 
the CR density around young sources. Responsible for 
these anisotropics are turbulent field modes with varia- 
tion scales much larger than the Larmor radius of CRs 
which mimic a regular field. This effect can be con- 
firmed via the observation of irregular gamma-ray emis- 
sions around CR sources. 

Cosmic rays in magnetic turbulence. — Since we want 
to test effects beyond the diffusion approximation, we 
propagate individual CRs in turbulent magnetic fields 
-B(fc) oc exp(— ikx) using the numerical code described 
in [lfj, |ll I . The validity of this code was checked re- 
producing earlier results from [l2|, EH • We assume that 
the spectrum V(k) of magnetic field fluctuations is static 
and follows a power-law, V(k) cx k~ a . The former as- 
sumption is justified, because the changes introduced 
by a finite Alfven velocity are small for the considered 
time scales. We fix the mean magnetic field strength 
£? 2 ms = (B 2 (x)) as B rnls = 4/xG, normalising B Tms for 
fluctuations bounded by l m - m = 1 AU and Z max = 150 pc. 
In the numerical simulations, we choose i m ; n sufficiently 
small compared to the CR Larmor radius Rr,. We also 
adopt an isotropic spectrum of fluctuations, as we want 
to demonstrate that even in this case CRs diffuse initially 
anisotropically. 

The spectral index of the turbulent GMF is only 
weakly constrained, and both Kolmogorov (a = 5/3) 
and Kraichnan (a = 3/2) spectra are consistent with 
observations 0, [3] ■ As our results do not vary much be- 
tween these two cases, we present them for a Kolmogorov 
spectrum. The turbulence is expected to have a Bohm 
spectrum (a = 1) only close to shocks, where efficient CR 
acceleration requires a diffusion coefficient D(E) ~ R L . 
As we are interested in time-scales when CRs have al- 
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FIG. 1: Eigenvalues di (solid lines) of the diffusion tensor 
Dij = (xiXj) / (2t) together with the average diffusion coeffi- 
cient D (dashed line) as function of time t. For -B rms = 4/iG, 
Z max = 150 pc, a = 5/3 and CR energy E = 10 15 eV. 



ready escaped from the acceleration zone and diffuse on 
scales 0(100 pc), we do not address here the question of 
CR diffusion in the shock region. This allows us also to 
neglect the backreaction of CRs on the turbulent field, 
discussed e.g. in [l5|, which modifies the magnetic field 
only in a thin layer in front of the shock. We also only 
consider CRs with E > 1 TeV, which are those relevant 
for gamma ray observations above 100 GcV. 

CR diffusion and diffusion tensor — Investigating the 
propagation of CRs close to their source requires to in- 
ject them localised in space, say at x = 0, and to 
propagate them in a given concrete realization of the 
turbulent field. We find that diffusion can be strongly 
anisotropic in a specific field realization, as long as the 



distance between CRs and the source is < L 



This 



anisotropy is washed out averaging over many realiza- 
tions of the turbulent field. Thus the correct procedure 
to calculate the diffusion tensor in this case is to compute 
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' for N particles (labeled by the 



(a) (a) 
ij _ 2Nt L-ia=\"i X j 

subscript a) and injected at x = in one single realiza- 
tion b. For each of the M realizations one diagonalizcs 
and finds its eigenvalues d\ b '. Then one averages 
the ordered eigenvalues, < d% < d$ , over the M 

realizations, di = -h YlbLi <H ■ 

In Fig. [1] we show the three eigenvalues di of 
as function of time for the case of CRs with energy 
E = 10 15 eV. We used M = 10 realizations, propagat- 
ing for each N = 10 4 particles. At early times t < t*, 
the diffusion tensor is strongly anisotropic, and the ra- 
tio ds I d\ between its largest and smallest eigenvalue can 
reach a factor of a few hundreds in some turbulent field 
realizations, while it is on average around a factor of a 
few tens. For t > t*, CRs propagate more and more 
isotropically, approaching the predictions of the diffusion 
approximation for a purely isotropic turbulent field. 

Modes with large scale variations contain most of 



the power in Kolmogorov turbulence, compared to the 
smaller scale variation modes which arc responsible for 
the diffusion of TeV-PeV CRs. Particles see modes with 
large variation scales as local uniform fields and diffuse 
therefore anisotropically. Once a sufficient fraction of 
particles moved beyond \x\ ~ l mra , the anisotropics of 
different "cells" of size ^ ax are averaged out and the CR 
densities around sources tend towards the limit predicted 
by the diffusion approximation. Adding a large scale reg- 
ular field Bq on top of the turbulence would increase the 
anisotropics. In this case, CRs arc known to diffuse faster 
in the direction of Bq than in the perpendicular direction. 

The dashed line in Fig. [T] represents the average diffu- 
sion coefficient D = D^, with defined as 
= l/(QNt) Yja=i x {a) -x^ for the magnetic field re- 
alization b. For the largest times numerically reachable, 
t = 10 5 yr, the eigenvalues in Fig. Q] and the average D 
approach a common value, D ss 5.5 x 10 28 cm 2 /s. Inter- 
estingly, 2 x 10 4 yr corresponds to (r 2 ) x l 2 = ^/6Dt ps ^ max 
valid in the isotropic limit. We verified that the limiting 
value for the average D agrees with the value of the diffu- 
sion coefficient computed using CRs with random start- 
ing positions. It is consistent, among others, with the 
computations of [l2| for pure random fields. The value 
of D(E) in our Galaxy is currently only known within a 
factor ps 50 at E ~ 10 GcV, and 5.5 x 10 28 cm 2 /s is in 
the acceptable range extrapolating D(E) to E = 10 15 eV 
for a Kolmogorov spectrum 1(| TtI- It is a factor ps 4 
smaller than the value used in Ref. [5j. 

Earlier works 0, El, IH HJ did not report anisotropic 
diffusion for several reasons: Diffusion coefficients were 
computed averaging over several configurations, with 
random initial positions for particles, or the considered 
space or time scales were too large. For instance, [3] 
and 12| find isotropic diffusion in the limit of a vanish- 
ing uniform field, Bq <g B rms . In both works the diffu- 
sion coefficient are calculated averaging over many reali- 
sations of the turbulent field. As a result, any anisotropy 
is averaged out, if the random field is isotropic. We 
also stress that t» is much larger than the transition 
time Tdiff ~ AD/c 2 from the ballistic to the diffusive 
regime [18j [. 

Cosmic ray intensity and extrapolation to low E. — In 
the middle row of Fig. [2l we show the projection of the 
number density of 1 PeV CRs on an arbitrarily chosen 
plane of size 600 pc x 400 pc containing the injection point 
x = in the center. We consider here one given realiza- 
tion of the turbulent field, out of the ten used for Fig. [T] 
The diffusion is confirmed to be strongly anisotropic at 
early times, 500 yr (left) and 2000 yr (middle panel), 
and even strongly filamentary. At 7000 yr (right panel), 
the distribution of CRs slowly tends towards the spher- 
ical limit expected for true isotropic diffusion. Most of 
the nine other configurations display similar filamentary 
structures at early times. For a few of them, no thin fil- 
aments are visible, but the CR distribution around the 
source is still asymmetric, showing wind-like structures, 
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FIG. 2: Relative cosmic ray densities around their source projected in the panel planes, for energies E = 100 TeV (upper row), 
IPeV (middle row), lOPeV (lower row) and times t = 500 yr (left column), 2kyr (middle column), 7kyr (right column). Same 
field realization in each panel. Each panel corresponds to a 600 pc x 400 pc field-of-view, with the source located in the center. 



also visible for some of the observed sources [7J. 

The upper and lower rows of Fig. [2] present results for 
E = 100 TeV and E = lOPeV, respectively. A compari- 
son of the three rows shows that the period of anisotropic 
diffusion lasts longer for CRs with lower energy. For in- 
stance, the panel with t = 2 kyr and E = 1 PeV is very 
similar to the one with t = 7 kyr and E — 100 TeV, which 
suggests that the expected scaling t oc 1/D(E) oc E" 1 ^ 
also holds in the case of anisotropic diffusion. To deter- 
mine this scaling law more quantitatively, we plot the 
values of df"\ and as function of time for 

given realizations b in Fig. [3] and look for times with 
similar d^ /d^: Numerically we find t* oc E~"< with 
7 = 0.25-0.5, i.e. a value of 7 consistent with the theo- 
retical expectation. 

Transition time. — For the parameters of Fig. [T] one 
estimates t* ~ 10 4 yr. This value is mostly deter- 
mined by Z m axj and we find that the naive expectation 
t* 'Lx holds in a first approximation. Reducing Z max 
by a factor 6, to 25 pc, the transition happens at 200- 
300 ~ 10 4 /6 2 yr in all ten tested configurations. There- 
fore, our numerical results suggest that the diffusion ap- 



proximation predictions become valid at 

U ~ 10 4 yr (? ma x/150pc)' 3 (E/PeVy (B rms /4 M G) 7 

(1) 

with (3 ~ 2 and 7 = 0.25-0.5 for Kolmogorov turbulence. 
B Ims ~ 4fj,G corresponds to the estimated local value 
around the Earth from unpolarized synchrotron radia- 
tion data [l9j]. However, the large uncertainties of the 
turbulent field parameters, in particular of l max , imply 
that i* may differ significantly from 10 4 yr for PeV CRs. 

The spectral index a does not have a strong impact 
on i*. In contrast, it influences the ratio d^/di at early 
times. For a Bohm spectrum, we find that the anisotropy 
a t t < t* is reduced. Indeed, more power is concen- 
trated in small scale modes than in a Kolmogorov spec- 
trum, resulting in a better isotropization of PeV CRs. 
After averaging over five realizations of the field, wc 
find a factor d^/di m 3-4 at early times. Out of these 
five configurations, only one is Kolmogorov-like with 
dtp /di ~ 10 and visible filaments. The others are just 
slightly anisotropic. 

Gamma-ray emission. — High energy protons can scat- 
ter on protons of the interstellar gas, producing sec- 
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FIG. 3: Eigenvalues of the diffusion tensor Dij as function of time t for energies E = 100 TeV (left), 1 PeV (mi 
lOPeV (right). Same field realization as in Fig. [2] Dashed lines for the average diffusion coefficient ~ D. 
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FIG. 4: Left panel: Relative cosmic ray densities along the 
lines-of-sight as seen from a specific observer located 500 pc 
away from the source. E = 1 PeV and t = 1 kyr, for a given 
magnetic field realization. Box size is 20° x 20°, with the 
source at (0,0). Right panel: Corresponding relative surface 
brightness in 7-rays with energies i? 7 > 100 GeV. 



of magnetic field fluctuations, I < i max - The propaga- 
tion of such CRs close to their sources has to be studied 
in single realizations of the turbulent field. We showed 
that CRs diffuse anisotropically at early times t < i», 
with from Eq. {l}, leading to a filamentary structure 
of the CR density around their sources. Turbulent field 
modes with variation scales much larger than the Lar- 
mor radius of CRs are responsible for this anisotropic 
diffusion regime. This effect can explain the observa- 
tions of irregular gamma-ray halos [6|-|8| around CR pro- 
ton and electron sources. If CRs propagate distances 
I ft imax, these anisotropics are averaged out. CR den- 
sities around sources become isotropic and tend towards 
those expected from the diffusion approximation. 

GG acknowledges support from the Research Council 
of Norway through an Yggdrasil grant. 



ondaries which in turn decay into photons. We sim- 
ulate cross sections and the final state of proton- 
proton interactions using QGSJET-II [2(|, while we use 
SIBYLL 2.1 [2l[ for the subsequent decays of unstable 
particles. The emission of secondaries is strongly forward 
beamed and CRs in filaments have an anisotropic distri- 
bution of momenta. We may expect that the emitted 7- 
rays are a good tracker of the underlying CR anisotropics. 

For simplicity, wc assume a uniform gas density around 
the source. Then we place an observer at 500 pc distance 
from the source and integrate the photon flux emitted 
along the line-of-sight towards the observer. We model 
the observer as a sphere of 5pc, which is the small- 
est size providing reasonable statistics and introducing 
only a small amount of artificial 'fuzziness'. The result- 
ing source image is shown in the right panel of Fig. |4] 
The comparison to the corresponding CR intensity (left 
panel) shows that the latter can be used to predict the 
shape of the gamma ray halo. Note that gamma rays 
emitted via Compton scattering by electrons would also 
display such anisotropic patterns. 

Conclusions. — We studied the diffusion of TeV-PeV 
CRs on scales I smaller or comparable to the largest scales 
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